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Abstract. We study the statistical properties of the time-averaged mean-square 
displacements (TAMSD). This is a standard non-local quadratic functional for inferring 
the diffusion coefficient from an individual random trajectory of a diffusing tracer in 
single-particle tracking experiments. For Brownian motion, we derive an exact formula 
for the Laplace transform of the probability density of the TAMSD by mapping the 
original problem onto chains of coupled harmonic oscillators. From this formula, we 
deduce the first four cumulant moments of the TAMSD, the asymptotic behavior 
of the probability density and its accurate approximation by a generalized Gamma 
distribution. 
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1. Introduction 

The statistical properties of Brownian motion play a crucial role in disciplines as various 
as mathematics, physics, chemistry, biology, engineering, economics and ecology [THS]. 
In many cases, an observation or an experiment can be repeated a large number of 
times to get a representative statistics of outcomes. In other situations, one measures 
an observable which reflects a cumulative effect of a large number of molecules (e.g., the 
macroscopic nuclear magnetic resonance signal formed by an extremely large number of 
nuclei [5]). In both cases, a system can be characterized by ensemble-averaged quantities. 
For instance, the ensemble average (or expectation) of the mean-square displacement 
(MSD) of a freely diffusing particle, (X 2 (t)), yields a measure of the diffusion coefficient 
D through the Einstein's relation: (X 2 (t)) = 2Dt. However, there are other situations 
when few or even single stochastic trajectory X(t) is available that makes the ensemble 
averaging inaccurate or impossible. The common examples are stock prices in finance [6] 
or trajectories of macromolecules inside living cells acquired by single particle tracking 
(SPT) techniques [3UTHT6]. In the former case, the observation (i.e., a price time series) 
is evidently unique for each asset. In the latter case, although an experiment can be 
repeated, the ensemble average may still be problematic, especially for tracers that move 
in spatially heterogeneous and time-evaluating media such as living cells. In order to 
infer the dynamical information from such individual trajectories, one needs therefore 
to replace ensemble averages by time averages. For instance, the diffusion coefficient is 
often inferred by considering the time-averaged MSD (TAMSD) at a lag time t over a 
sample trajectory of duration T: 



Although this is still a random variable, the time average reduces its fluctuations around 
the mean which, for Brownian motion, is (xt,r) — 2Dt. 

The statistical properties of quadratic functionals of a Gaussian process have been 
intensively studied. In mathematical statistics, several series representations of the 
probability density of a general quadratic form of a Gaussian process have been proposed, 
e.g. a mixture of chi-squared distributions pTHiM] . a power series expansion [251126] . 
and a series expansion over Laguerre polynomials [2TH29] (see also [301131] for a review). 
For biophysical applications, Qian and co-workers analyzed the TAMSD of a discrete 
random walk which is known to have a Gamma (or chi-squared) distribution only at 
the smallest lag time [7]. The distribution of diffusion coefficients was also studied via 
Monte Carlo simulations by Saxton 0|9] (see also [101132] for experimental data). More 
recently, Boyer and Dean investigated the distribution of several estimators of diffusion 
coefficients for Brownian motion [33J. They reduced the analysis of least-squares and 
maximum likelihood estimates to the distribution of a local quadratic functional of 
Brownian motion. Mapping this problem onto an imaginary time Schrodinger equation, 
they derived explicit analytical results for the distribution of the diffusion coefficient 
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estimator. In a separate paper, Boyer and co-workers discussed the optimal way to 
extract diffusion coefficients from single trajectories of Brownian motion and showed the 
superior efficiency of the maximum likelihood estimator over least-squares estimates [34J . 
An optimal Bayesian method for quantifying biomolecule diffusivity was presented by 
Voisinne et al. [35J. 

A different approach to the analysis of inferring schemes with quadratic functionals 
was developed by Grebenkov [361137] . First, the explicit formulas for the mean and 
variance of the TAMSD (and another quadratic functional) were derived for a large class 
of Gaussian processes governed by a generalized Langevin dynamics [36]. The formula 
for the variance allows one to estimate the spread of measurements and to choose an 
appropriate sample duration T. Second, a matrix representation for the characteristic 
function of Xt,T was used in order to analyze the probability density p(z) of the TAMSD 
of a general discrete Gaussian process [H7J- Moreover, an empirical approximation for 
p(z) by a generalized Gamma distribution (GGD) was proposed: 

/ x zV ~ X ( a z\ 
V{z) = ; f=- exp , 2 

with three parameters a > 0, b > and !/£R, and K v (z) is the modified Bessel function 
of the second kind (in the limit a — > 0, one retrieves the Gamma distribution ^-fj^p—)- 
Although the high accuracy of the GGD approximation (J2]) was confirmed numerically 
for a discrete Brownian motion, the theoretical status of this approximation remained 
unclear. In addition, the dependence of the parameters a, b and v on the lag time t was 
not analyzed. 

In this paper, we rigorously characterize the TAMSD of Brownian motion. In 
Sect. El we introduce a functional integral representation for a Laplace transform <p(s) 
of the probability density p(z) which is a cornerstone for all consecutive derivations. 
For t > T/2, it leads to a simple closed formula for ip(s). For t < T/2, an exact 
representation of (f(s) involves a determinant of an explicit matrix of size [T/t] x [T/t], 
[T/t] being an integer part of T/t. In Sect. [3j we discuss some practical consequences of 
these theoretical results. In particular, we derive the explicit formulas for the first four 
cumulant moments of the TAMSD, as well as the asymptotic behavior of the probability 
density p(z) at small and large z. For instance, we show that p(z) decays as e~ a l z at small 
z that justifies the use of a GGD as a convenient approximation for p(z). It also explains 
a remarkable accuracy of this approximation which correctly reproduces the asymptotic 
behavior of p(z) at both small and large z. Finally, we derive the asymptotic relations 
for the parameters a, b and v of the GGD as functions of the lag time t by matching 
the first moments of the true and approximate distributions. These relations allow one 
to get an accurate approximation of the probability density p(z) for a given lag time t 
that opens a number of practical applications for single-particle tracking experiments. 
These findings and future perspectives are summarized in Conclusion. 
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2. Theoretical results 

Throughout this section, we consider X(u) to be Brownian motion with mean zero 
and covariance (X(u),X(v)) = mm{u,v}, with diffusion coefficient D set to 1/2. The 
sample duration T is fixed to 1, while the lag time t e [0, 1] is dimensionless. We are 
interested in the probability density pt{z) = (5(z — xt)) of the random variable Xt which 
denotes the TAMSD ofX{u): 

l-t 

1 



o 



We aim at calculating the Laplace transform <ft(s) of the probability density Pt{z) 

oo 

cpt(s) = (exp(-sxt)) = / dze- sz Pt (z), (4) 



o 



Once tpt(s) is obtained, one can retrieve the probability density Pt(z) either through 
the inverse Laplace transform of ift{s), or through the inverse Fourier transform of the 
characteristic function <f>t(k) = (e lkxt ) = (p t (—ik) (although both ways are formally 
equivalent, Fourier transforms are more convenient for numerical implementation). 
From a theoretical point of view, the knowledge of <ft( s ) is therefore fully equivalent 
to that of Pt(z) (see Sect. [3] for practical consequences). 

We start by expressing the Laplace transform in Eq. (j3J) as a functional integral 
over trajectories X(u) of a Brownian motion pfldT] 

X(l)=xi i 1 l-t \ 

Ms) = jdx± J VX{u) exp -i J (d u X(u)) 2 - ^ J (X(u + t)- X{u)f , (5) 

R X(Q)=0 V / 

where we assume that Brownian motion starts from the origin: X(0) = (anyway, the 
TAMSD is translation invariant and thus the above functional integral is independent of 
the starting point). The next step would be to use the Feynman-Kac formula [4|[T8] and 
map the functional integral to a propagator of a certain quantum mechanical problem, 
to compute the propagator and to evaluate ^Pt(s). This method was successfully used 
by Boyer and Dean for local quadratic functionals of Brownian motion [33] • However, 
in our problem different parts of the trajectory, at times u and u + t, are coupled, 
and the action is non-local, forbiding a direct use of the Feynman-Kac formalism. To 
overcome this difficulty, we partition the trajectory into pieces and relabel them as new 
Brownian walkers in order to reduce it to an effectively many-body but local problem. 
The latter is not surprising since a system with memory is often equivalent to some 
many-body system, for example dynamics of a classical spin-glass or quantum mean- 
field spin-glass are typically mapped onto single particle problems with memory [T9l 20j. 
The transformation appearing in our problem is similar in spirit, but uses the inverse 
mapping to turn a non-local single particle problem into a local many-body problem. 
We illustrate this in detail for the case t > \ where the mapping and the computation 
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Figure 1. (Color online) Splitting of a trajectory into blocks, (a) For t > i, 
there are only three blocks: two coupled parts Yq and Y\ and one free part Zq; (b) 
When t = (m = 2,3,4,...), a trajectory is split into m equal blocks yielding a 
chain of m successively coupled oscillators Yq, F m _i; (c) When — ^ < t < 
(for some m = 2,3,4,...), the time interval [0,1] is covered by m subintervals of 
duration t, plus a smaller interval [1 — mt, 1] of duration 5 = 1 — mt. In order to split 
the interval [0, 1] regularly, each subinterval is divided into two parts, of durations 5 
and t — 6, respectively. As a consequence, the trajectory is mapped onto two chains 
of coupled oscillators {Yfc} and {Zk}- The two chains {Yfe} and {Zk} are coupled 
through boundary conditions for successive parts of the trajectory: Y k {5) — Z k (0) 7 
k = 0, ...,m — 1 which express continuity of the original trajectory X(u). 



are relatively simple. The case t < | is a straightforward generalization which is treated 
in the next Subsection. 



2.1. Case t > \ 

The partitioning of the trajectory in the case t < ~ is illustrated on Fig. [TJl. We split 
the whole trajectory X{u) into three parts: 

Y (u)=X(u) uE[0,l-i], 
Z (u) =X(u + l-t) ue[0,2t-l], 
Yi(u) = X(u + t) u G [0,1 -t]. 

The original trajectory is then replaced by a collection of two harmonically coupled 
particles and a free walker. Boundary conditions express continuity of the original 
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y (l-t)=a:i-t Z (2t-l)=x t Yi(l-t)=si 

<Pt(s) = J dx 1 _ t dx t dx 1 J VY J VZ J VY 1 

R 3 Vo(0)=0 Z o (0)=xi_ t Yi(0)=xt 

2t-l l-t 1-i 

,2 



X 



exp [~\jdu (d u Z ) 2 ~\jdu [(d u Y ) 2 + (fl^) 2 ] - |- | d« (y x - y ) S 





Note that the subscripts of x±- t , x t and Xi are just notations for these integration 
variables. We have defined for convenience 

2 2s 

9 = —f 

The original non-locality is translated into harmonic coupling of particles Y , Yi, the 
overall action is now local. Since there are just two particles involved, this case is easy 
to resolve: Y , Y\ are decoupled in a standard way by defining their normal modes A , 
A, 

A (u) = ^=(n(«) + Y 1 (u)), Ax(«) = -^(Yo(u) - Yi(u)), 

where Aq represents the center of mass motion and A\ describes harmonic oscillations 
around the center of mass. Boundary conditions for Aq and A\ are: 

-L(y (0) + y 1 (0)) = -^(0 + x t ), 

) = ^o(l ~t) + Y 1 (l- t)) = -L(xi_ t + xi), 

-L(y (o)-y 1 (o)) = -^(o-x t ), 
) = -L(y (i - 1) - - 1)) = ^(xi_ t - m). 

Therefore we get a system of two free particles Zq, Aq and one harmonic oscillator A\\ 

- J 4i(l-t)=ai,i_ t / i /.l-t 9 /-l-t 



°0,0 


= A)(o) 


a o,i-t 


= - 


°1,0 


= A(o) 


ai,i-t 


= ^i(i- 



<p t (s) = / dx 1 _ t dx t dx 1 / VA l exp -- / du(9 tt Ai) 2 - ^- / duA? x 

J 7Ai(0)=oi o V ^ ./0 -Wo / 

R 3 

pZ(2t-i)=x t / y r 2t ~ l \ fAo(l-t)=ao,i- t / -y [l-t \ 

/ VZ exp[-- du(d u Z ) 2 ) VAoexpl-- du(d u A ) 2 ) . 

JZ(0)=x 1 -t \ Z JO J JA (0)=ao,o V 2 ^0 / 

The Feynman-Kac formula can now be applied to each particle separately expressing 
<ft{s) as a product of propagators of respective quantum problems: 

<ft(s) = J dx 1 _ t dx t dx 1 G g (a 1A - t , l-t \ a lfi , 0) G (a ,i-t, 1 - £ | a ,o, 0) G (x t , t \ x^ t , 1 



i.e. two free propagators G and one propagator G g for a harmonic oscillator: 

2 / 

- pvn I — 

vM 1 - 1 2 ) 



G 9 (x,t\y,u) = - 7 ==== T e W [--(x + y) - -(x - y) ) , (6) 
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where q = e u ^ 9 and a — (1 — q)/(l + q). In the limit g — > 0, one retrieves the diffusive 
propagator of a free walker: 

Go (x, t | !/ , M ) = -^i= e X p(-|^g). (7) 
Pluging expressions for propagators, one gets 

(p t (s) = [ dxi-tdxtdxi = =- exp (-^f(a 10 + a 1A _ t ) 2 - -j-(a 10 - a x ,i- t ) 2 ) x 
J - Q ) V 4 4a y 

/ (a ,i-t-ao,o) 2 \ / (xt-xi-tf 

I 2 (l-t) J ^ I 2(2t-l) 
X v . — ^- X 



^2vr(l - t) v/27r(2t - 1) 

where 1 — q)/(l + q). Substituting the expressions for ao,o, ao,i-*, 

01,0, ai,i-i and evaluating directly the Gaussian integrals, we get the following closed 
formula 

-1/2 




^ = U+VaC 1 -*)^— T TT=^ • (8) 

This formula is one of the main results of the paper. 

2.2. Case t < \ 

Although the analysis for the case t < \ relies on the same ideas and it is a 
straightforward generalization of the above calculus, the computation is much more 
involved because the number of particles is growing linearly to infinity as t — > 0. 
Technically one has to distinguish subcases of < t < ^ for all integer m > 1. 
This is a reflection of the fact that the smaller t is, the more images t + ku with 
integer k can be placed inside the interval [0,1], yielding X(u) coupled to X(u + t), 
which is coupled to X(u + 2t) and so forth. As previously, the original problem with 
self-interaction maps onto an exactly solvable model of coupled harmonic oscillators. 

For discrete lag times t = — with m = 2,3,4, the whole time interval [0, 1] is 
covered by m equal subintervals so that the trajectory can be mapped onto a single 
chain of m coupled harmonic oscillators as shown on Fig. [Tb. Using this mapping, we 



derive in Appendix A.l the exact formula for the Laplace transform <pt{s): 



'm-\ ^ \ZA~~ \ ^ 1 



where = 2(1 — cos(— )), with k = l,...,m — 1. The matrix A of size m x m, 
which represents the coupling of harmonic oscillators, is given by Eqs. (1A. 12j) or (1A. 13j) . 
Although Eq. ()9]) provides an exact solution for the original problem, it has a rather 
formal nature, as evaluation of det(A) for large matrix sizes m (i.e. small times t) is 
quite complicated since A is a generic matrix with no particular simplifying structure 
(see 



Appendix A. 3 ). This is the main difficulty encountered in the computation of <ft( s ) 



for small values of t. 
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In the generic case < t < ^ (for some m = 2,3,4,...), the interval [0,1] 
is covered by m subintervals of duration t, plus a smaller interval [mt, 1] of duration 
5 = 1— mt < t. In order to split the interval [0, 1] regularly, each subinterval is divided 
into two parts, of durations 5 and t — 5. This covering maps the whole trajectory onto 
two interacting chains of m + 1 and m harmonic oscillators as illustrated on Fig. [lb. In 



Appendix A. 2 we derive the following exact formula 



, , l t~t Sgy/xje \ ' fyr (t-5)gy/Xj^ \ ' 1 . . 

rl " l ii sinh(^VS) J \l\smh((t-5)gVX'k)J y/6^(t - 5)™ det(A) ' 



where = 2(1 — cos(^j)), with k — 1, m. The matrix A of size (2m +1) x (2m + 1), 
which represents the coupling of harmonic oscillators, is given by Eq. ( 1A.17I) . Once 
again, this solution is exact though still quite formal, as the dependence of <ft(s) on s 
and t is partly "hidden" in the determinant of the matrix A. 

3. Discussion 

Quite surprisingly, a theoretical analysis of the statistical properties of the TAMSD 
turned out to be a challenging problem, even in the simplest case of Brownian 
motion. This reflects the non-locality of the involved quadratic functional and the 
many-body character of the underlying quantum system. It is also worth stressing 
that a "resolution" of this problem would mean different things for theoreticians and 
experimentalists. In the previous Section, we have fully resolved the problem from a 
theoretical point of view, by providing the exact formulas ©, ©, (fTOj) for the Laplace 
transform <^t{s) of the probability density. However, this solution remains somewhat 
formal since an explicit "shape" of the probability density is still unknown, even for the 
simplest case t > |. In this Section, we discuss some practical consequences of these 
results. 

3.1. Moments 

3.1.1. Case t > ~. Using the exact solution (jEJ), one can compute the moments of the 
TAMSD: 

"* = (- 1 )*(^( a ))^' ( U ) 
from which the cumulant moments are 

«i = t, (12) 

K2 = ^(m 2 -6t + i), (13) 

k 3 = — (286t 3 -216t 2 + 54t-4), (14) 
15 
2 

k 4 = (8065i 4 - 8220t 3 + 3206t 2 - 572t + 41). (15) 

105 
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Figure 2. (Color online) (Left) The reduced variance K2/K1, skewness kz/k 2 and 
kurtosis k^/k^ of the TAMSD \t as functions of the lag time t. Symbols present the 
values obtained from an exact matrix formula for a discrete random walk (see Eq. (6) 
from [37], the sample length N = 1000), while lines refer to the analytical formulas 
(Ull), O, 0U for t > 1/2 and OH), {33]), JE) (or (2D)) for * < 1/2. (Right) The 
relative error between the analytical formulas and their counterparts for a discrete 
random walk. The vertical dotted lines stand at t = 1/3 and t = 1/2. All the formulas 
are accurate for the whole region (an increase of the relative error at small t is related 
to intrinsic differences between discrete and continuous processes). 



Here, k% is the mean value, k 2 the variance, while the skewness and kurtosis are given 
as Ko,j K<j 2 and k±/ k\, respectively. We note that the expression for the variance k 2 
reproduces the result by Qian et al. which was obtained for a discrete random walk [7] . 
In turn, the formulas for the third and fourth moments are derived for the first time. 



3.1.2. Case t < 1/2. For the case 1/3 < t < 1/2 (m = 2), the explicit formula for 



tpt(s) is provided in Appendix A. 3 Taking the derivatives with respect to s, we obtain 
the first cumulant moments: 



K 3 = — 



t, 

1 t 3 (4-5t) 
3 (1-t) 2 ' 

2 t 5 (33 - 47t) 
15 (1-t) 3 

2 



(16) 
(17) 



K4 



105 



30913t 8 - 86272t 7 + 102060t° 

3 



(19) 



68040t 5 + 28350t 4 - 7560t 3 + 1260t 2 - 120t + 5 



As the explicit formula ( 1T0|) for ipt{s) is too cumbersome for t < 1/3, another 
representation ( IB. 21) for the cumulant moments was deduced in Appendix B In 
particular, we checked that Eqs. (fTTj) and (jTBl) for k 2 and K3 are also valid for t < 1/3. 



(20) 



In turn, the formula for k 4 for t < 1/3 becomes 

8t 7 



k 4 



105(1 - 1) 4 



(302 -473t) (t < 1/3). 
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Figure |2] shows the reduced variance K2/K1, skewness k^/k^ 2 and kurtosis k^/k 2 , of 
the TAMSD as functions of the lag time t. The values obtained from an exact matrix 
formula for a discrete random walk (see Eq. (6) from |37J, the sample length N = 1000), 
are compared to the analytical formulas ffT3l). (PT4D, (fTol) for t > 1/2 and (TTT1). ffT8l). (fT9l) 
(or (120]) ) for t < 1/2. Notice an excellent agreement between the curves. 

3.2. Small z asymptotic behavior 

3.2.1. Case t > | In the limit of large s, Eq. ([8]) can be approximated as 



2e~v /s ( 1 -*) 
<p t (s) » — = (s>l). 

Using the asymptotic relation for the inverse Laplace transform L -1 , 

K+1/2 

L -l[ s » e -2V55] ~ ^^^-3/2^/^ (21) 
'7T 



with k = —1/4 and a = (1 — 1)/4, one derives the asymptotic behavior of the probability 
density pt{z) at small z: 

9 e -(l-t)/(4z) 

p t (z)n = (z^0). (22) 



l-t 



In other words, the asymptotic decay of ift{ s ) as e~ 2 ^ at large s implies the sharp 
decay of pt{z) as e~ a ^ z at small z. This is the first rigorous justification for using a GGD 
in Eq. (j2J) to approximate the probability density Pt(z)- 

3.2.2. Case t<| In this case, we have to analyze the asymptotic behavior of Eqs. 
(Q or (fTUl) . For large s, one can replace sinh(,z) by e z /2 (with exponentially small 
corrections) to get in both cases 



exp ( -y/s- 





' J2(l-t) , 

<*(*) « V 7= 7 (« » !)• (23) 



where the functions r](t) and /t(g) are explicitly given in Appendix A The asymptotic 



behavior of (pt(s) at large s is mainly determined by the exponential term, while h{g) 



turns out to be an algebraic correction. In Appendix C we derive a representation of 



It(g) in the form h{g) = c + cig, where c and C\ become independent of g for large g 
(or s). As a consequence, the asymptotic behavior of ipt(s) reads as 



exp ( -y/s- 



' J2{l-t) 

Neglecting c and using again Eq. (l2Tj) . one deduces the asymptotic behavior of pt(z): 



Time-averaged MSD of Brownian motion 



11 



in which Ci(t) is given by Eq. (1(140 . Similar to the case t > |, the probability density 
sharply decays (as e~ a ^ z ) when z — > 0. 

3.3. Large z asymptotic behavior 

Since the inverse Laplace transform can be written as the Bromwich integral in the 
complex plane, the singularities of <ft{s) for s G C play an important role. In 
particular, the singularity with the largest (negative) real part, which is called the 
abscissa of convergence of the Laplace transform, determines the asymptotic behavior 
of the probability density pt(z) at large z [38J. 

For t > |, one can easily check that all singularities of (pt(s) are located on the 
negative real axis. In fact, taking s = — a 2 /(l — t), one can rewrite Eq. (jSJ) as 

<p t (s) = — (l -ft" 1 atana)~ 1/2 , (25) 
cos a v 

where /3 = (1— i) / (3i— 1) is a positive parameter between and 1 (given that | < t < 1). 
This function has two sets of singularities at points and ol\~ which are solutions of 
the equations 

cosafc = 0, a^tanafc = (3 (k — 1, 2, 3, ...). 

The smallest solution a>i of the second equation can be determined perturbatively in 
powers of /3: 

«? = /3-^ 2 + ^/3 3 + 0(/3 4 ). 
The related singularity is thus located at 

„ = -jL = _ ( i _ Jf±^ + Jil^iL + o ( (i - «)■; 



1-t V. 3 ^- 1 3(3t-l) 2 45(3t-l) 3 
Note that ot\ varies from at t — 1 to 0.8603... at t — 0.5 so that \si\ varies from 0.5 
at t — 1 to 1.4803 at t — 0.5. At the same time, the smallest solution a\ = tt/2 is 
independent of t and is always larger than a±. 
In close vicinity of s\, one has 



from which the inverse Laplace transform yields an asymptotic approximation 



2(/3 2 + a 2 ) exp(-z/b) , 

"■-^^ (g+A^-D ^r 1 (2>>1) - (26) 



with 6 = l/|si| = (1 -t)/a\. 

The analysis can be extended to the case t < |. For instance, when t — 1/m, one 
can substitute s = —(1 — t)a 2 /(2t) into Eq. (jUJ) in order to locate the singularities of 
(pt{s) either as solutions of sin(afc-y/Afc) = (with k — 1, ...,m — 1), or as solutions 
«fc of det(A) = 0. For a general situation ^-j- < t < ^ for some m = 2,3,4, the 
singularties of <ft{s) can be determined in a similar way. In both cases, the large z 
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Figure 3. (Color online) The probability density pt(z) (solid line), its small z 
asymptotic relation ([2H]) (shown by circles), and large z asymptotic relation (|!26| (shown 
by crosses), at lag times t = 0.9, t = 0.75 and t — 0.5, in semilog (left) and loglog 
(right) scales. 



asymptotic behavior is determined by the largest negative root Si (i.e., the root with 
the smallest absolute value |si|) of the equation det(A) = 0. Given the square root type 
of singularity in Eq. (Tl0|) . the exponential asymptotic decay of pt(z) at large z in the 
form similar to Eq. ( 126]) is expected to be valid for any lag time t. 

The accuracy of the asymptotic Eq. (122]) for small z and Eq. ( 126]) for large z is 
illustrated on Fig. [3] For all lag times, Eqs. (|22]) .(l2^j) correctly describe the asymptotic 
behavior of pt{z) at small z. Moreover, for larger t (e.g., t = 0.9 and t = 0.75), Eq. ([22]) 
turns out to reproduce the whole shape of pt(z), except for large z. The deviations at 
large z are expected because the exponential decay of pt(z) is not captured in the small 
z asymptotic relation ( 122]) . 
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3.4- Limiting case t = 1 

In the limit t — > 1, one gets ^fi(s) = 1/V 1 + 2s, from which 

= (27) 

V27T2; 

We retrieved therefore the limiting behavior of the discrete case for the lag time 
n = N — 1 [7J[37]. Since the density p\{z) describes the distribution of a squared 
Gaussian variable Xi = -^ 2 (1)> fi{ s ) an d then Eq. (1271) could be directly obtained by 
integrating the Gaussian probability density -^ e ~ x2 ^ 2 f° r -^"(1) with e~ sx2 . 

The function pi(z) diverges at z = and monotonously decreases on the positive 
semi-axis z > 0, in sharp contract to pt{z) (with t < 1), which must have a maximum, 
as being equal to at z = and z — > oo. In the limit t — >■ 1, the most probable value of 
Xt (i-e., the position of the maximum) approaches 0, while the mean value approaches 
1. This also illustrates the fact that p±(z) has the largest skewness among p t (z), as seen 
on Fig. H 

3.5. Limiting case t — > 

In the opposite limit of t — > it is harder to get analytic behavior from theoretical 
results, as the limit describes the situation where the number of particles m grows 
to infinity, which, combined with complicated boundary conditions for normal modes, 
makes a direct evaluation unfeasible. We provide here a few qualitative arguments, 



while a more rigorous analysis is reported in Appendix D Since the mean and variance 
of Xt vanish as t — > 0, it is convenient to rescale Xt by t, in order to fix the mean value 
to 1. For small t, the rescaled TAMSD can be approximated as 



~ = *i ~ _ YV 2 

t m ' 

k=l 



where m = [1/t] and (k are independent Gaussian displacements (standing for 
X(to + t) — X(to), with to = kt) with mean zero and variance 1. The above sum is 
known to have a chi-squared (or Gamma) distribution with m degrees of freedom |39j, 
from which 



<M S J 



1 + 2s/m)- m/2 ~ (1 + 2st)- 1/{2t) . (28) 



In the limit t — )■ 0, we get <^o( s ) = e ~ s an d Po(z) = 5(z — 1) as expected since Xt/t 
concentrates around 1 in this limit. Although the limiting results for (f>o{s) and po(z) 
are correct, an approximation ( 128]) of (^t(s) for small t by a Gamma distribution is not 
accurate. 

3.6. Dimensional units 



All the above results were obtained for D = 1/2 and T = 1, with t being a dimensionless 
lag time between and 1. Dimensional units can be easily incorporated by replacing s 
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by 2DTs, z by z/(2DT) and t by t/T in the formulas for (pt(s) and pt{z). For instance, 
Eq. ((SJ) reads for a given diffusion coefficient -D and sample duration T as 

2e-v^M^) / ^ 3t _ T i_ e -2v^(7-7)N ~ 

<A (s) = . 1 + J2Ds{T - t) — 

In the formulas for the cumulant moments, one has to replace t by t/T and multiply 
the right-hand side by (2DT) n , with n being the order of the moment. 



3.7. Discrete versus continuous cases 

As we discussed in Introduction, the TAMSD is commonly employed for inferring the 
diffusion coefficient from an individual random trajectory of a diffusing tracer. For this 
purpose, one records positions of a tracer at successive times, with a fixed time step r. 
An individual trajectory of the tracer is therefore discretized, and the TAMSD becomes 

^ N-n 

Xn,N = ~rf ( X k+n ~ X kf (29) 

k=l 

where = X(kr) is the position of the tracer at time kr, n — t/r is the discrete lag 
time, and iV = T/r is the number of points in the sample. In the limit r — > 0, the 
discrete TAMSD is expected to converge to the continuous one defined by Eq. (0Q). The 
probability distribution of the TAMSD for a discrete Gaussian process was investigated 
in [37]. For a given covariance matrix C (and zero mean), one can use the classical 
matrix representation 

<fn,N( S ) = (exp(-SXn,v)) = -==!==, (30) 

^/det(I + sMC) 

where the matrix M represents the discrete quadratic form in Eq. (1291) . When 
the sample length iV is smaller than few thousands, the Schur decomposition of the 
matrix MC allows one to rapidly compute the Laplace transform ip n ^(s) or the 
characteristic function 4> n ,N(k) = <p n ,N(—ik), from which the probability density p(z) 
can be reconstructed through the inverse Fourier transform [37J. 

Although Eq. (I3"0~|) provides an exact solution for v ? n,A r (s) and an efficient numerical 
way for computing p(z) for moderate N, the dependence on the parameters (e.g., the lag 
time) remains elusive. In contrast, Eq. OH]) for Brownian motion is explicit and allows 
one to derive the asymptotic behavior and the dependence on the lag time. Since Xn,N 
approaches \t,r in the limit r — > 0, Eq. (EJ) is expected to be an accurate approximation 
to v 9 n,v(s) for large enough N. This point is illustrated on Fig. H] in which (ft(s) 
for t = 0.5 is compared to (Pn/2,n(s), for different sample lengths N. The curves for 
N = 8, 16, 32 exhibit deviations at large s, while the curves for larger iV = 64, 128 are 
barely distinguishable from ipt(s) at the present range of s values. As a consequence, one 
can use the explicit formula (JHJ) for discrete samples with iV > 50. At the same time, the 
functions y n ,Ar(s) and (pt(s) will always be different for large enough s. In fact, for a fixed 
N, v ? n,v('5) exhibits a power law asymptotic decay (as det(I + sMC) is a polynomial 
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Figure 4. (Color online) (Left) Comparison between ipt(s) for Brownian motion 
(with t = 0.5) from Eq. (J5J (thick black line) and fNt,N(s) for discrete random walks 
with different sample length N = 8,16,32,64,128, with n = tN. The curves for 
N = 64, 128 are barely distinguishable from the theoretical one. (Right) The relative 
error between (pt{s) and </Jjvt,Jv(s). As expected, the relative error increases with s and 
it is larger for smaller N. 



in s), while <y3 f (s) decays exponentially fast. However, the distinction between these 
behaviors may appear for so large s, for which both functions become already extremely 
small so that a distinction would be irrelevant for practical applications. In summary, 
Eq. OH]) and its extensions for t < 1/2 can in practice be used for studying the TAMSD 
of discretely sampled positions of a freely diffusing particle. 



3.8. Relation to generalized Gamma distribution 

As shown numerically in [37], the probability density Pt(z) for a discrete random walk 
can be accurately approximated by a generalized Gamma distribution (GGD) from Eq. 
(EJ). Its Laplace transform reads as 



K v {2^Jb) V ; 



On one hand, its comparison to Eqs. (JBJ), fl9]), (TTOj) immediately shows that a GGD 
is not an exact distribution for \t (except for the trivial case t = 1). On the other 
hand, we have seen earlier that the GGD correctly reproduces the asymptotic behavior 
at small and large z, and it turns out to be a remarkably accurate approximation 
of pt{z) (see Fig. [6]). Moreover, a simple form of Eq. ([2]) is particularly attractive 
for practical purposes (e.g., inferring statistics of diffusion coefficients by maximum 
likelihood methods, estimating the confidence intervals, etc.). For practical uses of this 
approximation, it is important to relate the parameters a, b and v to the lag time t. 



In Appendix E we derive the asymptotic relations for the parameters a, b and v in 



the limit t — > by matching the first three moments: 



3V1T59 ! , 

xv ~ — r 1 + c w 0.638 r 1 - 0.537, 32 

160 

1 63 

- ~ — r 1 + u w 0.394 - 0.714, (33) 
£ 160 




Figure 5. (Color online) The parameters a, b and v (and c = 2^/a/b) of the GGD 
as functions of the lag time t. Full circles are obtained by a numerical resolution of 
Eqs. (|E.2j) . (|E.3[) and solid lines represent the short-time asymptotic formulas (|32[) . 

(22), (2D, (ESD for t < i. 



+ a t ^0.178 - 0.105*, (34) 




l + v/TT^ 2 320 

— t 2 + 6 3 t 3 « 1.749 t 2 + 2.3t 3 , (35) 
183 

/ Tl59 



where c = and x = ^ ~ 1.6211. Here, the leading terms were derived 

analytically, while the values of the next-order corrections u , Co, a and b were obtained 
by fitting the numerical curves of these quantities. Although they could in principle be 
found from the higher-order terms, such analysis goes beyond the scope of the paper. 
The asymptotic relations ( 1321) . ( 1331) . ( 1341) . ( 1351) turn out to be accurate for the whole 
range of t < ~, as shown on Fig. [5j 

We check the accuracy of the GGD with the parameters a, b and v from the 
short-time asymptotic formulas (133]) . (|34|) . (|35|) by comparing this approximation to 
the probability density obtained numerically for a discrete random walk (Sect. 13. 71) . 
The latter is referred here as "exact" solution, although it was obtained by a numerical 
computation of the inverse Fourier transform of the exact characteristic function from 
Eq. ( 1301) . as explained in detail in [37]. In order to plot several curves onto the same 
figure, the TAMSD was rescaled by t to set the mean value to 1 for all lag times t. Figure 
|6]4eft) shows an excellent agreement between exact solutions and approximations by the 
GGD, for several lag times t from to ~. 

For larger lag times, the probability density Pt(z) becomes more and more skewed, 
with the maximum of Pt(z) getting much smaller than the mean value. As a consequence, 
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Figure 6. (Color online) The probability density Pt(z) of the TAMSD rescaled by t 
for several lag times t, in normal, semilog and loglog scales. Lines show exact solutions 
for a discrete random walk, while circles represent approximations by GGDs. For t < i 
(left panel), the parameters a, b and v are obtained from the asymptotic formulas (|33|) . 
(f34|) . (|35[) . Even for the limiting value t = \, the approximation remains accurate. For 
t > i (right panel), the parameters a, b and ^ are extracted from the best fit of Pt(z) 
by a GGD. 



positive-order moments are determined by the values of Pt(z) far from the maximum. 
Given the objective to approximate the function pt(z) around the maximum, one cannot 
therefore rely on matching positive-order moments, as done in Appendix E One needs 
therefore to identify other characteristics of the true and approximate distributions to 
be matched. One can see that even with the exact solution at hand, the analysis is 
intricate. 

At the same time, GGD remains an accurate approximation of the probability 
density pt{z) even for t > |, as illustrated on Fig. [6] (right panel). On the right 
panel, the parameters a, b and v were obtained by fitting the exact density pt{z) 
(obtained numerically) to a GGD. As a perspective, one may attempt to derive analytical 
asymptotic relations for these parameters for t > |. Note also that, for t > |, the small 
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z asymptotic formula (122]) is also an accurate approximation of Pt{z), as discussed in 
Sect. IO 

3.9. Two-dimensional and higher- dimensional case 

The above analysis was focused on a one- dimensional Brownian motion. In turn, many 
SPT techniques allow one to record the displacements of a tracer in two or three 
dimensions. Under the assumption of isotropic motion, one often considers the combined 
time-averaged MSD 

T-t d 

xl d T = j^~ t [ dt J2(Mt + t)- X t (t )) 2 , (36) 
o 1=1 

where Xi(t) is the tracer's position at time t along the ith coordinate. For d- dimensional 
Brownian motion, each process Xi(t) is a one-dimensional Brownian motion independent 
from the others so that the Laplace transform <Pt(s) for the random variable xfr * s 
simply [<^ t (s)] d . One can then easily extend the above ID results to higher-dimensional 
cases. In particular, the probability density p(z) would exhibit similar asymptotic 
behavior at small and large z, with deviations coming from power law corrections. 

Conclusion 

We have studied the statistical properties of the time-averaged mean-square 
displacements of Brownian motion. Using the standard tools of functional integrals, 
the original problem was mapped onto an exactly solvable model of coupled harmonic 
oscillators. When the lag time t was larger than a half of the trajectory length T 
(t > T/2), a simple closed formula (jHJ) for the Laplace transform (pt(s) of the probability 
density pt(z) was derived. In the opposite case t < T/2, the analysis was quite similar, 
although the resulting exact formulas were less explicit. This reflects the many-body 
character of the original problem. In the special case t = 1/m with m = 2,3,4,..., one 
dealt with a chain of coupled harmonic oscillators and derived Eq. (Q for <ft{s)- The 
general situation corresponded to two chains of coupled harmonic oscillators, for which 
Eq. ffTUl) for <pt(s) was deduced. From these formulas, the probability density p t (z) 
can be formally reconstructed either through the inverse Laplace transform of (pt(s), or 
through the inverse Fourier transform of the characteristic function <f>t(k) = ift(-ik) (the 
second option being more convenient for numerical implementation). To our knowledge, 
the probability distribution of the TAMSD of Brownian motion was derived for the first 
time. 

From the exact formulas for (pt{s), we deduced the first four cumulant moments 
of the TAMSD. While the first and second moments (mean and variance) were 
reported earlier, the third and fourth moments (related to skewness and kurtosis of the 
distribution) have been obtained for the first time. We also analyzed the asymptotic 
behavior of the probability density Pt{z) at small and large arguments. In particular, 



Time-averaged MSD of Brownian motion 19 

we have shown the sharp decay of pt(z) as e~ a / z at small z. This behavior justified 
the use of a generalized Gamma distribution as an approximation for pt{z) suggested 
in [37]. The fact that the GGD correctly reproduces the asymptotic behavior of the 
probability density at both small and large z explained a remarkable accuracy of this 
approximation, illustrated by Fig. El Finally, we expressed the parameters a, b and 
v of a GGD as functions of the lag time t by matching the first moments of the true 
and approximate distributions. The resulting asymptotic relations were shown to be 
accurate for t < T/2 so that one could easily construct an accurate approximation for 
the probability density Pt(z) of the TAMSD for a given lag time. For larger t, the 
distribution was too skewed so that the first moments did not really control the shape 
of Pt(z) around the maximum. In that case, one had to resort to a simple approximation 
of pt(z) provided by the asymptotic form (l22~j) (i.e., by a GGD with a = (1 — £)/4, b = 
and v = 0.5). However, as illustrated on Fig. El better approximations by a GGD could 
be constructed. Their explicit derivation requires further analysis. 

From a practical point of view, an accurate approximation of the probability 
density Pt(z) by a simple function (such as a generalized Gamma distribution) opens 
several interesting perspectives for a more reliable analysis of random trajectories of 
diffusing tracers acquired in single-particle tracking experiments. In particular, one 
can develop efficient inferring schemes through maximum likelihood estimation, or 
accurately quantify the confidence intervals of measured diffusion coefficients. While 
the present work was focused on Brownian motion, future investigations may challenge 
the statistical properties of the TAMSD for more sophisticated processes (e.g., restricted 
diffusion, fractional Brownian motion or generalized Langevin equation with memory) 
which are known to be representative models for diffusive dynamics in viscoelastic media 
(e.g., living cells) or for time series in finance. 



Appendix A. Computation for t < | 

The computation for the case t < ~ is a straightforward generalization of the calculus 
from Sect. 12.11 For convenience, we consider separately the special case t = — in 
Appendix A.l and the generic case ^—^ < t < — in Appendix A. 2 



Appendix A.l. Special cases t = — 

For t — — with m = 2,3,4,..., the whole time interval [0,1] is covered by m equal 
subintervals so that the trajectory can be mapped onto a chain of m coupled harmonic 
oscillators as illustrated on Fig. [lb. We define oscillators Y k as 

Y k (u) = X(kt + u) ue{0,t), (k = 0...m - 1), 

so that the non-local part of the action in Eq. (JSJ becomes: 

l-t (k+l)t t 

du{X{u + t)- X{u)f = J2 M X ( U + t)~ x (u)f = J2 MYk+i(u) - Y k (u)) 2 , 

k=0 kt k=0 
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m - 1 '^(t)=JHl 



dy x ...dy m [ j / W fc exp 

k=0 JY k (p)=y k 



(d u Y k f - 9 - 



(Y k+1 - Y k f 



with the starting point y — (and y m = xi). The non-locality is therefore eliminated 
at the expense of many-body character of the resulting problem. To advance further we 
need to diagonalize the chain of coupled oscillators with open boundary conditions: 



TO— 1 „t TO — 2 TO — 1 „ 

E = £ / (W 2 + 5 2 £ / (Y k+1 -Y k f = Y: / 

!,_n ^0 u- n JO u _ n JO 



m—l 



{duYkY + g^YkLkM 



1=0 



k=0 " u k=0 " u k=0 

where the matrix L is a discrete Laplacian of size m x m with von Neumann boundary 
conditions: 

/ 1 -1 ••• \ 



-12-10 
0-12-1 
0-12 
-1 













\ o o o o ••• -l i / 

The normal modes which diagonalize S are defined as 



TO — 1 



A k = J2^kiYi (k = 0,...,m-l), 



1=0 



where is an orthogonal matrix with columns formed by eigenvectors of L: L-^ = A. 
It is straightfoward to compute the eigenvalues and eigenvectors of this matrix: 

7T/V 



Afci = 5ki\ k , A fc = 2 1 - cos 



m 



(k = 0...m - 1), 



• nk(l+l) ■ nkl 

sin — - — - — sin — 

* W = ; m g - (k = 1...771 - 1, Z = 0...771 - 1), 



0/ 



^ m (l_cos^) 

(Z = 0...771 - 1). 



(A.l) 
(A.2) 



We obtain the diagonalized form as 

TO-l pt 



s = E / k^) 2 + 



(A.3) 



(note that for k = the second term vanishes as A = 0), with the boundary conditions 
for A k translated from the boundary conditions for Y k as 

to — 1 m—l 

a k0 = A fc (0) = E *fci^(0) = E ( k = 0-m - 1), 



TO — 1 



i = l 

m—l 



Ofct 



= A k (t) = *kiYi(t) = *kiyi+i (k = 0..m - 1). 



(A.4) 



1=0 
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These linear relations convert the original coordinates yi, . . . , y m into new variables a^o 
and akt- In the first relation, the sum runs from I = 1 as yo = and the columns of the 
matrix ^ are shifted to the left, the last column being filled with zeros (as there is no 
dependence on y m ). This shift can be represented by a shift matrix S (of size m x m) 
as * ■ S: 



/ 





. 


. 


\ 




1 


. 


. 










1 . 


. 





\ 





. 


. 1 


o / 



Using the shift matrix, we can write 

m— 1 

ciko = ■ S]fc;j/;+i (k = 0...m - 1). 

1=0 



(A.5) 



Once we have the diagonalized form flA.31) . the Feynman-Kac formula is applicable 
and the Laplace transform (pt(s) is again given as a product of m — 1 propagators of 
harmonic oscillators and one propagator of a free particle (corresponding to Ao = 0): 

/m— 1 
dyi ■ ■ ■ dy m G (a t,t\aoo,0) G gy ^(a kt , t\a k0 , 0). 



k=l 



We plug in the propagators as given by Eq. (jSJ) and find the following expression: 



<M S J 



dyi--- dy r , 



exp 



x — exp 



(a>ot ~ a oo) 
2t 

4 

1 - gfc 

1 + Qk 

The formula for <£>t(s) can be rewritten as 

exp (-\gr]{t)) 1 



(A.6) 



O'kt + O'kO) : 

4ai 



where 



9* = e 



-tgVXk 



Ctk 



m— 1 

n (i + 

fc=i 



(A.7) 



(A.8) 



in which 



m— 1 



r/(t) =t^ v / ^ = t(ctan(£-)-l 



(A.9) 



fe=i 
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where y = (yi, y m ) T , and A = V T • UV is the matrix of size m x m which combines 
both the structure of the quadratic form in Eq. flA.6j) and the linear relations flA.4j) . 
( 1A.5I) . represented respectively by matrices U and V of sizes 2m x 2m and 2m x m: 



U 



u+ 


u 


u 


u+ 



V=(-^-l, (A.11) 



where U are diagonal matrices of size m x m with the elements 

u t = ^9V^k(a k + 1/afc) = ta J^^p^ (k = l...m-l), u t = \' 

u k = lgVh(®k ~ V«fc) = . r7~T^ (A; = l...m - 1), u = — -. 

The elements Uq can be seen as the limit of the generic formulas with Ao — > 0. 
Performing the matrix multiplication yields 

A = V T UV = * T U + * + S T * T U * + * T U *S + S T * T U + *S. (A.12) 

The diagonal matrices U contain the eigenvalues X k of the original matrix L and ^ is 
formed by its eigenvectors, implying a simple matrix representation for A: 



A = - S T 9VZ ^ - dVZ ^ S + S T 9VZ ^ S. (A.13) 

tanh(ptvL) smh(gtvL) sinh(gi(;vL) tanh(gtvL) 

This expression is convenient for a theoretical analysis, while the original representation 
(IA.12I) is better suited for rapid numerical construction of A. 

The computation of the Gaussian integral in Eq. (1A.10I) yields 

1 = gg) ^ ffr 1 V2 

V 7 ^ Vdet(A) vij 

This relation, together with Eq. (IA.8j) . completes the computation of ift{s) for special 
cases t = 1/m (with m = 2,3,4,...) and provides a solution to the problem. This 
solution can alternatively be written in the form 0. 



Appendix A. 2. Generic situation 

For the generic case < t < — (for some m = 2, 3, 4, ...), the interval [0, 1] is covered 
by m subintervals of duration t, plus a smaller interval [mt, 1] of duration 5=1 —mt < t. 
In order to split the interval [0, 1] regularly, one divides each subinterval in two parts, 
of durations 5 and t — 5, respectively: 

[0, 1] = [0, 5) U [6, t] U • • - U [tm-utm-l [tm-1 + S, t m ) U [t m , 1], 

where t k = kt. This covering maps the whole trajectory onto two interacting chains of 
m + 1 and m harmonic oscillators as illustrated on Fig. [Tb: 

Y k {u) = X{t k + u) ue[0,S\, (k = 0...m), 
Z k (u) = X(t k + u) ue[6,t], (fc = 0...m-l). 
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Here m — [l/t], the integer part of 1/1 One gets therefore 

r m ~ 1 rY k (5)=z k 

(p t (s) = / dy 1 ...dy m dz ...dz m TT / VY k exp 

fc=0 JY k (0)=y k 



1 g 2 * 



{d u Y k Y-»- I (Y k+1 -Y k ) 2 



'm-2 ,Z k (t-5)=y k+1 

■ < n / vZk exp 

"Y m (S)=Zm 

x / exp 

'Y m (0)=y m 



^ (d u Z k ) 2 - y / (Z,,+ | - Z/,.j 



2J 



l(*)=»ti 



PZ m _i exp 



m-l(&)=Z m -l 



(with 2/o — 0)- There is no direct interaction between the two chains, but the boundary 
conditions for the chains are mixed since Y k (5) = Z k (0) and Z k (t — 5) = Y k+1 (0) to 
ensure the continuity of the original Brownian trajectory. Introducing normal modes: 

m m—l 

A k = ®kiYi, B k = J2 *kiZi, 

1=0 1=0 

this problem transforms into a set of uncoupled oscillators. The \I/ is an orthogonal 
matrix formed by the eigenvectors of the discrete Laplacian L of size (m + 1) x (m + 1) 
(tilde sign refers to quantities related to the matrix L of size (m + 1) x (m + 1), in order 
to distinguish them from similar quantities corresponding to L of size m x m). The 
boundary conditions read for each k — . . . m 



k=i 

m—l 
1=0 

We obtain then 



a k s = 2j ^klZl 
1=0 

m 

=1 



i-iVi- 



(A.15) 
(A.16) 



/( m ^ ( m—l 

dyi...dy m dz ...dz m < JJ G^^-(a kS , 5\a k0 , 0) I <^ JJ G g ^{b kt , t\b kS , 6) 

B2i7i+l Kk=0 ) \ k=0 



Note that two factors with k = correspond to free diffusion propagators (as A = A 
0). Substituting Eq. (JHJ) for the propagators, one gets Gaussian integrals: 

1 / (ao<5-aoo) 2 (bot-bos) 2 " 



<Pt{s) 



dy 1 ...dy m dz ...dz 



exp 



2m+l 



X 11 ^d-g) CXP 



27Tv/5(t - 5) 



2(t - 5) 



fc=i 

m—l 



i-r (^VAD 1/2 e- ( *- 5)s ^ /2 
x exp 



4 
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. \2 9v^k f n2 



4at 



(6« + 6. 



(frfct — fefed)' 
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The argument of the exponential function can be represented as a quadratic form 



i(a T Ua): 
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b m -i,t 
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Q>m6 
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where a is the (4m + 2)-dimensional vector, U is a (4m + 2) x (4m + 2) matrix with the 
elements: 



■it. 



-g\/X k (a k ±l/a k ) (k = l...m), 



it, 



g-\/\ k (a k ± l/a fc ) (/c = l...m - 1), w 



±1 

T' 

± 
o 



±1 



The intermediate variables aoo, ^oar ■ ■ forming the vector a, are expressed through Eqs. 
(IA.15p . (IA.16P which can be written in a matrix form a = Vx, where x is (2m + 1)- 
dimensional vector x = y 2 , ■■■,y m , z o, z i, z m ) T , and matrix V has a block structure 



«m0 




^01 • 


^0m 











*11 • 


* lm 
























*00 • 


^0,m-l 











*m-l,0 • 


^m— l,m— 1 



















^O.TO-l 










^10 • 


■ *l,m-l 


*lm 








*m0 • 


^ m,m— 1 


^ mm 








*00 • 


^0,m-l 











^m-1,0 • 


• ^m-l.m-l 






/ Vi \ 

Urn 
ZO 



We obtain therefore the quadratic form |(x T V T UVx) = |(x T Ax), where 

A = V T UV 
is the (2m + 1) x (2m + 1) matrix. 



(A.17) 
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exact 

approximate 




Figure Al. (Color online) The function rj(t) and its approximation (IA.22|) for 
t < 0.5. For t > 0.5, one has the exact relation -q{t) — -\/2(l — t). 



The computation of Gaussian integrals yields an exact representation: 



where 



<Pt[s) 
Qk 

v(t) 

1 



vfe=l 



1 



+ qk 



Sgy/X. 



-(t-S)gVT k 



'm-l 

n 

,k=l 



Ok 



exp(-f^)) 
1 + ^ J ^/ijg) 



(A.18) 



Oik 



1 - gfc 
1 + qk 

l - 

1 + qk 



At = 2 ( 1 — cos 



m + 1 



7rA; 

A*. = 2 ( 1 — cos ( — 

m 



m— 1 



k=l 



k=l 



'n#)"Trif >w 



One can alternatively rewrite flA.18j) in the form ( ITU]) . 

Using explicit expressions for A&, A& and formula (1A.9j) . we get 



(A.19) 
(A.20) 



7T 



n(t) = ( ctan . 



— ctan ( | 

V4m/ 



(A.21) 



£ ( 1 + m ctan 



7T 



, (m + l)ctan , . 
4(m + 1) / V4m/ 



If £ is small, then we can approximate m ~ from which the series expansion of 
ctan(x) leads to 



r l (t)^--t-^-t 2 + 0(t i ) (t < 0.5). 

7T 1Z 



(A.22) 



This expansion accurately approximates the function rj(t) for t < 0.5, as illustrated on 
Fig. |ATJ For t > 0.5, one has m = 1, and Eq. (IA.21I) simplifies to r)(t) = y/2(l - t), 
while Eq. (jA.18|l is reduced to Eq. 
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Appendix A. 3. The case 1/3 < t < 1/2 

For the case 1/3 < t < 1/2 (m = 2), one performs the computation analytically and 
gets the following result 

_ exp(-f [(l-2t)(l + v^) + (3t-l)v^]) 

(l + ftXi + feXi + fcVcFF^ ' { } 

where 

Coffi") = 7=77 f 9 + 4v / 3<Si<52 + l0V6a 2 cni + a\a\ + 16a|a? + 6v / 2aiai«2 
576 V 

ci((?) = ^6«i + — j^t\a\a\ + 6ai«2 Q; i + 2v / 3a2Cti + ~^ ai ~^ 6v / 6ai«2cn + 2v / 3«i<52 

It is worth recalling that the "coefficients" c and c 1; as well as gi, g 2 , 9i> «i, ct2, «i ah 
depend on g which enters through an exponential function. As t gets smaller, expressions 
are getting even more cumbersome. 

Appendix A. 4- Exact formulas for t = 1/3 and t = 1/4 



For the case t = 1/m with small m = 2,3,4, the matrix A from Eq. (1A.13|) has a 



small size so that one can compute det(A) explicitly. For instance, one has for m = 3, 

9 2 
\2a\a 2 



g2 |- 

det(A) = — 4gv / 3«i + Aga 2 a\ + a\a\\[?> + 9^ + \2ct\a 2 



and for m = 4 

g 3 r 

det(A) = — — Aga 2 a\oi\\Plv\ + \2ga\azot\v\\[2 — 32a 2 asUiV2 

64a!a 2 «3 L 

+ 128aia 2 v / 2^i + 12v / 25'^a 3 + 2Sga 1 V2v l + %go? 2 a\a. x v x 
+ lQga\azCt\ui + 160aia2^i + 96z/ia2«3 + 40gaiz/i — 16gz/ia 3 
— 16a + 16aialV2a 3 + \1 ga 2 a\\Pl 

+ 17gaja 2 V2 + 64 + 2^a\a\ + 2Aa 2 2 a\ + 2Aga\a 2 - 2Aga 2 a\ + I2ga 2 a 3 a x 



where v x = ^/Jh/2 = sin(vr/8) = V 2 - v^/2. For lar ger m, explicit formulas are much 
more cumbersome. 



Appendix B. Cumulant moments 

The simple explicit formula ([8]) for <^t{s) allows one to easily compute high-order 
cumulant moments of Xt for t > |. The analysis of the opposite case t < | is more 
difficult, especially for small t. For this reason, we present an alternative approach to 
evaluation of these moments. 

We start by establishing a formal representation of (pt{s) as a determinant of some 
integral operator. Using the standard Gaussian identity 

e ~X 2 /2 _ f e -y2/ 2+ixy 

J V2^ 

R 
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one can rewrite Eq. §5§ by introducing an auxiliary Gaussian field Y: 

(l-t \ X(l)=ici 

-\ f Y 2 (u) j 'dx! J VX{u) 
/ R X(0)=0 

l-t 



xexp f-i J( du X(u)y + i(J^ J 



1/2 

Y(u)(X(u + t)-X(u)) 



Splitting the trajectory X(u) into a linear part tx\ and Brownian bridge X(u) started 
from and arrived at at time t — 1, one integrates over the whole path X(u) to get 

l-t l-t l-t 



¥>t{s) 



J VY(u)exp j ~ J duY 2 {u) - J du J dv Y{u)G{u,v)Y(v) 





where G(u, v) is a combination of the two-point correlators C(u, v) of Brownian bridge 
X{u) 

C(u, v) — (1 — v)u — (u — v)Q(u — v), 

G(u, v) = t 2 + C{u + t,v + t) - C{u + t,v)- C{u, v + t) + C(u, v), 

where Q(x) is the Heaviside function: Q(x) = 1 for x > and otherwise. One gets 
then 

G(u, v) = (u-v + t)Q(u -v + t) + (u-v- t)Q(u -v-t)-2(u- v)Q(u - v), (B.l) 

which also reduces to t — \u — v\ for \u — v\ < t, and otherwise. 

Evaluating Gaussian integrals over the auxiliary field Y{u) yields 



det (1 + ^G) 



where G is an integral operator with the kernel G(u,v) acting from L2QO, 1 — t\) to 
L 2 ([0, 1 — t)). This is a direct extension of the classical representation (as Eq. ( 1301) ) 
for discrete Gaussian processes to the continuous case (see [37J). Rewriting the above 
formula as 

_ ln „ (s) = ^ ln ( 1 + ^)^f;I(^i)" w 

n=l 



one gets a representation of cumulant moments: 
More explicitly, one has 

i-t i-t 
2 n ~ l (n — 1)! f f 
f^n 7i Ivn / dui... I du n G Ul U2 G U2 ,v,3 ■•■ G Un lll , (B.2) 



with an explicit function G uv from ( IB. II) . 
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Appendix C.l. Case t = 1/m: representation of the determinant 



In this section, we only consider the case t = 1/m, with m = 2,3,4,..., although the 
general situation could be analyzed in a similar way. The matrix U from Eq. ( 1A.11I) 
can be split into two parts: U = Uo + gUi, where the matrix Uo has only 4 nonzero 
elements which do not depend on g: 



( 



Ur 



Un 



Ur 





ut 



\ o ... o o 

Given that — an d u o = ~' 1 



\ 






U 1 = p" 1 (U-U ) 



/ o 



An 



V T U V 



1/t = m, one finds 
\ 



V 









1 



/ 



One has therefore 

rn 

det(A) = det(A + g Ai) = det(Ai) det(^I + A^A ) = det(Ai) + a k ), 

k=l 

where A x = V T UiV and are the eigenvalues of the matrix Aj" 1 A . Since the matrix 
A has only one nonzero element, the product A^ 1 A is a matrix in which only the 
last column is nonzero. As a consequence, all the eigenvalues of this product are zero, 
except one, which is equal to (A]" 1 ) m _i im _i. We obtain therefore 

det(A) = det(A 1 )( 7 m - 1 (^ + (A^ 1 ) m _i )m _i) = g m ~\a + a l9 ), (C.l) 

where the coefficients an = det(Ai)(A^ 1 ) m _i i?Ti _i and a± = det(Ai) depend on g only 
through the coefficients in the matrix Aj. Denoting 

m— 1 

(I, i r a k ^ 



m2 



m—l 



n 

k=l 



0,1), 



one gets It(g) — cq + gc\ so that 

exp 



m— 1 

fc=i 



(1 + gi)...(l + q m -i)y/c + gci 
This representation is exact but quite formal as the dependence on g and t is still 
"hidden" in coefficients Cq and C\. 
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Figure CI. (Color online) The exact coefficient c\(t) for m = 2, 101 computed 
through Eq. ([(I2j) and its fit as 2£ 2 (0.213 - 0.145i 2 ). 



Appendix C.2. Case t = 1/m: large g asymptotic behavior 
For large g, all are close to 1, while are negligible so that 

C --vS^T (0»M = O,1), (C.2) 
where we used the identity 

m— 1 m— 1 , 

n v 7 ^ = n 2 sin ^) = ^ ( c - 3 ) 

fc=i fc=i 

Since the coefficients a; depend on g only through they rapidly converge to the 
limiting values. In the limit g — > oo, the matrix Ui gets a diagonal form 

from which 

Ai = V T UiV ~ V T v / AV + S T V T VAVS = VZ + S T VhS {g > 1). 

The limiting matrix on the right-hand side, which is explicit and independent of g, 
determines the limiting values ao and a±. Moreover, the explicit formulas for and ^ki 
allow one to get the elements of the matrix y/h (k,l = 0, ...,m — 1): 



Va. 











'hi 



m \cos(^) -cos(^M) cos (_^) _ cos (^±^ 



Using this representation, we obtain numerically an accurate approximation 
ai = det(Ai) ~ 2 m (0.213 - 0.145t 2 + ...), 

from which 

<*(*) - ~v|^I = 2t 2 (0.213 - 0.145t 2 + ...). (C.4) 
Its accuracy is illustrated on Fig. IC1I 
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Appendix C.3. General case 

In the general case ^jrj < t < ^ (for some m = 2,3, 4, ...), the analysis is quite similar 
but formulas are more cumbersome. The contribution h{g) from the Gaussian integral 
reads as 

h{9) = 6{t ~ P det(A). 

2 2m -V m "V m ( m + 1 ) 

By construction, the matrix A can be represented as A = Ao + gAi, where the matrix 

A is independent of g, while the matrix A x depends on g only through and As 

a consequence, Ai becomes independent of g for large g (or s). The determinant of the 

matrix Ao + <?Ai is in general a polynomial of g of degree 2m + 1. We conclude that the 

asymptotic behavior of <ft{s) in Eq. (fTOj) is essentially determined by the exponential 

function exp(— -y/i — ^= = ), with an algebraic prefactor from det(A). 

For the cases m = 1 and m = 2, we obtain that det(A) = g 2m ~ 1 {ao + a\g), with 
the coefficients ao and a\ depending on t. As a consequence, one gets h(g) = Co + cig, 
with the new coefficients Co and ci. We expect that this representation is valid for any 
m, but the related analysis is beyond the scope of the paper. 

Appendix D. Small g asymptotic behavior for t < \ 

We only consider the specific cases t = 1/m, with m = 2,3,4,.... One can split the 
explicit representation flA.13j) into two parts: A = H + Hi, with 

Ho = Hl = -S r gV V - + S T 9VZ ^ S. (D.l) 

tanh(ptvL) sinh(gtvL) sinh(gtvL) tanh(gtvL) 

One gets then det(A) = det(H ) det(I + H 1 H 1 ), from which the explicit computation 
of the determinant det(Ho) yields 

Up to this point, this is just another representation of the previous result. 

In the limit g — > 0, one can expand q\~ from Eq. (IA.7j) into Taylor series to get 

\fc=l / fe=l v 7 fe=l fe=l fc=l 

Given that 

m— 1 m— 1 

^A fe = 2(m-1), ^A 2 =6m-8, 
fc=i fe=i 



one finds 



rrl + cA V2 /Wv^ 1 s 2 t 3 (l-4t/3) 
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Note that the first term here exactly compensates the first factor in Eq. (ID.2[) . 

Now we turn to the analysis of the determinant det(I + Hq 1 H 1 ). Expanding 
expressions for Hq 1 and Hi from Eqs. ( ID. II) into Taylor series, one gets 

H x Hi = (S T S - S T - S) + ~ (2S T LS + S T L + SL - 2LS T S + 2LS T + 2Ls) 

-—( 2S T L 2 S + — (S T L 2 + L 2 S) + 10LS T LS + 5LS T L + 5LSl) + 0(z 3 ), 
90V 4 / 

where z = g 2 t 2 . Denoting J = I + S T S — S T — S, one has 



where 



Checking that det(J 
determinant, 



det(I + H^Hi) = det(J) det ^1 + zJ t + —J 2 + 0(z A ) 

Jx = - J" 1 (2S T LS + S T L + SL - 2LS T S + 2LS T + 2Ls) , 

J 2 = — J" 1 ( 2S r L 2 S + -(S T L 2 + L 2 S) + 10LS T LS + 5LS T L + 5LSl) . 

45 V 4 / 

= 1 and using a perturbative expansion for the second 



1, 



det (I + X) = 1 + tr(X) + ^[tr(X) 2 - tr(X 2 )] + ■ • • , 
one finds, up to z 2 : 

2 

det(I + H ^i) = 1 + ztr(Ji) + Z — (tr(J 2 ) + (tr(J x )) 2 - tr(J 2 )) + 0(z 3 ). 

Using the explicit form of the matrices Ji and J 2 , one can check that tr(Ji) = and 
compute exactly two other traces, from which 

^(2/3-t) 



detp + H^Hi) 



Bringing all these results together, we obtain 



0{s% 



<ft(s) 



cxpfj 2 t 3 (1 ~ 4</3) ' 



+ 0(s 3 ). 



(D.3) 



(D.4) 



One can see that the factor in front of e~ ts gives the second-order correction in the 
form s 2 t 3 . If one rescales the TAMSD Xt by t (to set the mean value to 1), the Laplace 
transform (fit(s') of the probability density for Xt/t is obtained by replacing st by s'. In 
the limit t — > 0, one can show that (f>t{s') converges to e~ s ', as expected. 



Appendix E. Comparison of moments 

Given that a generalized Gamma distribution turns out to be an accurate approximation 
for the probability density Pt{z), one can attempt to retrieve the parameters a, b and v 
of this distribution by matching the moments of the true and approximate distributions. 
For a GGD, one has 

Ha = (Xt) = c = 2v^76. (E.l) 
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Following Ref. [37J, one can solve the system of two nonlinear equations involving the 
first three moments: 

fi2 K v+2 {c)K v (c) /X3 K u+3 (c)K 2 (c) 

£ K 2 +1 (c) ' Kl +1 {c) ■ 1 • } 

Once the pair {c, v\ is found, one retrieves the other two parameters as 

6 = 20! *"®, v a = bc 2 /A (E.3) 
ciVn(c) 

(note that the first relation was misprinted in [37]). Using the explicit formulas (Tl2|) . 
(|T3|) . (JHJ) for the cumulant moments for t < |, one finds 



w 3 k 2 k 3 i(60 - 69t - 19t 2 ) 37, , 

q = l + 3^ + 4 = 1+ 1 — ^ ,. 3 ^ ^ 1 + 4t + -t 2 + 0(f). E.5 
/if k{ k\ 15(1 — t) 01 5 

We have solved numerically Eqs. (1E.2|) and found the dependence of the parameters a, 
b and v on the lag time t, as shown on Fig. O 

One can see that when t — > 0, both c and z/ increase inversely proportional to t. In 
this limit, one can approximate the modified Bessel function K v (yx) as 



K„(vx)~yfrfil 1 + Vl + x2 ) (l + x 2 )" 1 ^^. (E.6) 



,r 



Considering e — 1/V as a small parameter and assuming c/V to be a constant in the 
limit t — > 0, one takes c/V = x — + /3e 2 and gets 



K u+2 {c)K u (c) e f ax vl + 



2 



r- < - 1+ 7=T=f + 7TT 2A3/2 " 7TT 2^2 W + 0{e% (E.7) 

^ + i(c) vi + x 2 \{i + x 2 y/ 2 (l + x 2 ) 2 j 

K u+3 (c)K 2 (c) _ 3e f 3ax(2 + yTT^) 



Kl +l { c ) vTTx 2 V(l + x 2 )(l + v / T 

2(2 + x 2 ) v / rTx 2 + 3x 4 + Ax 2 + 4 



£ 2 + 0(e 3 ). 



(i + x 2 ) 2 (i + vTT^ 2 ) 2 

It is worth noting that these terms are independent of (3 (which appears in higher-order 
terms). 

Comparing these expansions to Eqs. (1E.4I) . (1E.5I) . one sees that e is of the order of 
t. Taking e = -ft + 5t 2 + 0(t 3 ), one immediately gets 7 = |a/1 + x 2 from the first-order 
terms. Substituting e = |a/1 + x 2 t + 5 t 2 into Eqs. ( IE. 71) . ( IE. 81) . one gets a system of 
two equations with three unknowns a, 5 and x: 

5 (ax vi + x 2 — 1\ 16 



vTT^ 2 Vv/TT^ 2 (1 + x 2 ) y 9 

35 /3ax(2 + vTTx 2 ) 2(2 + x 2 ) y/1 + x 2 + 3x 4 + 4x 2 + 4\ 16 _ 37 

VTTx 2 + V (l + vTT^ 1 ) 2 (i + x 2 )(i + vTT^ 2 ) 2 / "9" ~ ~5~' 

These equations are linear in a and 5. Expressing a from the first equation and 
substituting into the second one yields, after algebraic simplifications, a simple equation 
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on x which turns out to be independent of 5: 2x 2 — 38 + + x 2 (21x 2 — 38) = 0. Its 
solution is 

*= 2^ „ 1.8211. 
21 

As a consequence, one finds the leading terms for v and c and then for a and b as t — > 
which are summarized in Eqs. ( 1331) . ( l34l) . (J35]). 
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